setwd("C:/Users/geoff/OneDrive")

#Load packages
library("plyr")
library("dplyr")
library("rmarkdown")
library("randomizr")
library("multiwayvcov")
library("miceadds")

#Read in dataset
df <- read.csv("JEPS Final Data.csv")
summary(df)
str(df)

#Set block variable to class "factor"
df$block <- as.factor(df$block)

#VOTING HISTORY CATEGORIES
#Create subsets by voting history
df.highprop<-df[df$highprop==1,]
df.medhigh<-df[df$medhigh==1,]
df.medlow<-df[df$medlow==1,]
df.lowprop<-df[df$lowprop==1,]
df.newreg<-df[df$newreg==1,]
df.unreg<-df[df$unreg==1,]

#TABLE AB1. DESCRIPTION OF SAMPLE
#Calculate demographic statistics
table(df$racename)
12787/(8958+12787+14154+3411+1+2103) #30.88 percent Asian/Pacific Islander
8958/(8958+12787+14154+3411+1+2103) #21.63 percent Black
3411/(8958+12787+14154+3411+1+2103) #8.24 percent Latinx
14154/(8958+12787+14154+3411+1+2103) #34.18 percent White
table(df$age)
df$youth<-NA
df$youth[df$age<35]<-1
df$youth[df$age>=35]<-0
table(df$youth)
24718/(24718+16696) #59.69 percent under 35
table(df$unreg)
3674/(3674+37740) #8.87 percent
table(df$newreg)
2575/(2575+38839) #6.22 percent
table(df$lowprop)
8697/(8697+29438) #22.81 percent
table(df$medlow)
18495/(18495+19640) #48.50 percent
table(df$medhigh)
7391/(7391+30744) #19.38 percent
table(df$highprop)
3552/(3552+34583) #9.31 percent

#FIG. AH1. RANDOMIZATION DISTRIBUTION OF ITT, SHARP NULL
#Estimate Model 1 without block fixed effects to generate point estimate for randomization inference
model1.noblock <- miceadds::lm.cluster( data = df , formula = voucher ~ assignment, 
                                cluster = "householdID" )
summary(model1.noblock)

set.seed(101)
sims <- 10000
itt.store.sn <- rep(NA,sims)

# Set stratification, balance variables, and creative fields
strat_vars <- c("votepropensity","multivoterhousehold","wacanmember", "momcat", "reachable_cellcat", "reachable_emailcat")
bal_vars <- c("agescale","sex","racename","incomescore","educationlevel")
#Subset to required fields
household.sub <- subset(df, select = c(strat_vars, bal_vars,"householdID","randomizationrecords","voucher","block"))
#Remove all but household representatives from the subset
household.sub <- household.sub[household.sub$randomizationrecords==1,]

# Compile randomization distribution of the ITT under the sharp null
for (i in 1:sims){
  household.sub$assign_treat<-NA
  household.sub$assign_treat<-block_ra(blocks=household.sub$block,prob=0.5)
  bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
  p<-(coef(summary(bal)))[,4]
  p<-as.numeric(p)
  while(p[1]<0.1|p[2]<0.1|p[3]<0.1|p[4]<0.1|p[5]<0.1|p[6]<0.1|p[7]<0.1|p[8]<0.1|p[9]<0.1|
        p[10]<0.1|p[11]<0.1|p[12]<0.1|p[13]<0.1|p[14]<0.1|p[15]<0.1|p[16]<0.1|p[17]<0.1|
        p[18]<0.1|p[19]<0.1|p[20]<0.1|p[21]<0.1|p[22]<0.1|p[23]<0.1){
    household.sub$assign_treat<-NA
    household.sub$assign_treat<-block_ra(blocks=household.sub$block)
    bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
    p<-(coef(summary(bal)))[,4]
    p<-as.numeric(p)
  }
  
  ##Link subset back to main DF
  #Create assignment variable, assign duplicates to the values of their designated housemate
  df$D <- NA
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==1], 1, df$D)
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==0], 0, df$D)
  #Drop households with 8 or more voters
  df<-df[!df$householdID==24341,]
  df<-df[!df$householdID==27505,]
  df<-df[!df$householdID==1341,]
  df<-df[!df$householdID==58,]
  df<-df[!df$householdID==10117,]
  df<-df[!df$householdID==15890,]
  df<-df[!df$householdID==17965,]
  df<-df[!df$householdID==13328,]
  df<-df[!df$householdID==3609,]
  df<-df[!df$householdID==28910,]
  df<-df[!df$householdID==371,]
  df<-df[!df$householdID==1327,]
  df<-df[!df$householdID==13751,]
  df<-df[!df$householdID==27364,]
  df<-df[!df$householdID==4832,]
  df<-df[!df$householdID==13641,]
  df<-df[!df$householdID==25362,]
  df<-df[!df$householdID==9485,]
  df<-df[!df$householdID==10331,]
  df<-df[!df$householdID==18711,]
  df<-df[!df$householdID==19278,]
  #Estimate ITT under sharp null iteratively, storing each estimate
  itt.store.sn[i] <- mean(df$voucher[df$D==1],na.rm=TRUE) - 
    mean(df$voucher[df$D==0],na.rm=TRUE)
}
hist(itt.store.sn*100, main="Fig. AH1. Randomization Distribution of ITT, Sharp Null",breaks=50, xlab="ITT, percentage points")
abline(v=mean(0.003697596*100), col="red", lty=3)

sum((itt.store.sn)>(0.003697596))/sims


#FIG. AH2. RANDOMIZATION DISTRIBUTION OF ITT FROM INTERACTIVE MODEL, SHARP NULL

#Estimate Model 3 without block fixed effects to generate point estimate for randomization inference
model3.noblock <- miceadds::lm.cluster( data = df , formula = voucher ~ assignment+highprop+member+assignment*highprop+assignment*member, 
                                           cluster = "householdID" )
summary(model3.noblock)

#H2, H3: Repeat randomization inference for interactive LPM w/ categorical variable for high-propensity
# Compile randomization distribution of the ITT under the sharp null
set.seed(101)
sims <- 10000
itt.store.sn <- rep(NA,sims)

# Set stratification, balance variables, and creative fields
strat_vars <- c("votepropensity","multivoterhousehold","wacanmember", "momcat", "reachable_cellcat", "reachable_emailcat")
bal_vars <- c("agescale","sex","racename","incomescore","educationlevel")
#Subset to required fields
household.sub <- subset(df, select = c(strat_vars, bal_vars,"householdID","randomizationrecords","voucher","block"))
#Remove all but household representatives from the subset
household.sub <- household.sub[household.sub$randomizationrecords==1,]

for (i in 1:sims){
  household.sub$assign_treat<-NA
  household.sub$assign_treat<-block_ra(blocks=household.sub$block,prob=0.5)
  bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
  p<-(coef(summary(bal)))[,4]
  p<-as.numeric(p)
  while(p[1]<0.1|p[2]<0.1|p[3]<0.1|p[4]<0.1|p[5]<0.1|p[6]<0.1|p[7]<0.1|p[8]<0.1|p[9]<0.1|
        p[10]<0.1|p[11]<0.1|p[12]<0.1|p[13]<0.1|p[14]<0.1|p[15]<0.1|p[16]<0.1|p[17]<0.1|
        p[18]<0.1|p[19]<0.1|p[20]<0.1|p[21]<0.1|p[22]<0.1|p[23]<0.1){
    household.sub$assign_treat<-NA
    household.sub$assign_treat<-block_ra(blocks=household.sub$block)
    bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
    p<-(coef(summary(bal)))[,4]
    p<-as.numeric(p)
  }
  ##Link subset back to main DF
  #Create assignment variable, assign duplicates to the values of their designated housemate
  df$D <- NA
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==1], 1, df$D)
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==0], 0, df$D)
  #Drop households with 8 or more voters
  df<-df[!df$householdID==24341,]
  df<-df[!df$householdID==27505,]
  df<-df[!df$householdID==1341,]
  df<-df[!df$householdID==58,]
  df<-df[!df$householdID==10117,]
  df<-df[!df$householdID==15890,]
  df<-df[!df$householdID==17965,]
  df<-df[!df$householdID==13328,]
  df<-df[!df$householdID==3609,]
  df<-df[!df$householdID==28910,]
  df<-df[!df$householdID==371,]
  df<-df[!df$householdID==1327,]
  df<-df[!df$householdID==13751,]
  df<-df[!df$householdID==27364,]
  df<-df[!df$householdID==4832,]
  df<-df[!df$householdID==13641,]
  df<-df[!df$householdID==25362,]
  df<-df[!df$householdID==9485,]
  df<-df[!df$householdID==10331,]
  df<-df[!df$householdID==18711,]
  df<-df[!df$householdID==19278,]
  #Estimate ITT under sharp null iteratively, storing each estimate
  itt.store.sn[i] <- coef(lm(voucher~D+member+highprop+D*member+D*highprop, data=df, na.action=na.exclude))[6]
}
hist(itt.store.sn, main="Fig. AH2. Randomization Distribution of ITT from Interactive Model, Sharp Null",breaks=50, xlab="ITT, percentage points")
abline(v=0.011381598, col="red", lty=3)
mean(itt.store.sn)

sum((itt.store.sn)>0.011381598)/sims


#FIG. AH3. RAND. DIST. OF ITT ON HIGH-PROPENSITY VOTERS, SHARP NULL
#Estimate Model 8 without block fixed effects to generate point estimate for randomization inference
model8.noblock <- miceadds::lm.cluster( data = df.highprop , formula = voucher ~ assignment, 
                                        cluster = "householdID" )
summary(model8.noblock)

set.seed(101)
sims <- 10000
itt.store.sn <- rep(NA,sims)

# Set stratification, balance variables, and creative fields
strat_vars <- c("votepropensity","multivoterhousehold","wacanmember", "momcat", "reachable_cellcat", "reachable_emailcat")
bal_vars <- c("agescale","sex","racename","incomescore","educationlevel")
#Subset to required fields
household.sub <- subset(df, select = c(strat_vars, bal_vars,"householdID","randomizationrecords","voucher","block"))
#Remove all but household representatives from the subset
household.sub <- household.sub[household.sub$randomizationrecords==1,]

# Compile randomization distribution of the ITT under the sharp null
for (i in 1:sims){
  household.sub$assign_treat<-NA
  household.sub$assign_treat<-block_ra(blocks=household.sub$block,prob=0.5)
  bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
  p<-(coef(summary(bal)))[,4]
  p<-as.numeric(p)
  while(p[1]<0.1|p[2]<0.1|p[3]<0.1|p[4]<0.1|p[5]<0.1|p[6]<0.1|p[7]<0.1|p[8]<0.1|p[9]<0.1|
        p[10]<0.1|p[11]<0.1|p[12]<0.1|p[13]<0.1|p[14]<0.1|p[15]<0.1|p[16]<0.1|p[17]<0.1|
        p[18]<0.1|p[19]<0.1|p[20]<0.1|p[21]<0.1|p[22]<0.1|p[23]<0.1){
    household.sub$assign_treat<-NA
    household.sub$assign_treat<-block_ra(blocks=household.sub$block)
    bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
    p<-(coef(summary(bal)))[,4]
    p<-as.numeric(p)
  }
  
  ##Link subset back to main DF
  #Create assignment variable, assign duplicates to the values of their designated housemate
  df$D <- NA
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==1], 1, df$D)
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==0], 0, df$D)
  #Drop households with 8 or more voters
  df<-df[!df$householdID==24341,]
  df<-df[!df$householdID==27505,]
  df<-df[!df$householdID==1341,]
  df<-df[!df$householdID==58,]
  df<-df[!df$householdID==10117,]
  df<-df[!df$householdID==15890,]
  df<-df[!df$householdID==17965,]
  df<-df[!df$householdID==13328,]
  df<-df[!df$householdID==3609,]
  df<-df[!df$householdID==28910,]
  df<-df[!df$householdID==371,]
  df<-df[!df$householdID==1327,]
  df<-df[!df$householdID==13751,]
  df<-df[!df$householdID==27364,]
  df<-df[!df$householdID==4832,]
  df<-df[!df$householdID==13641,]
  df<-df[!df$householdID==25362,]
  df<-df[!df$householdID==9485,]
  df<-df[!df$householdID==10331,]
  df<-df[!df$householdID==18711,]
  df<-df[!df$householdID==19278,]
  #Estimate ITT under sharp null iteratively, storing each estimate
  df.highprop<-df[df$highprop==1,]
  itt.store.sn[i] <- mean(df.highprop$voucher[df.highprop$D==1],na.rm=TRUE) - 
    mean(df.highprop$voucher[df.highprop$D==0],na.rm=TRUE)
}
hist(itt.store.sn*100, main="Fig. AH3. Rand. Dist. of ITT on High-Propensity Voters, Sharp Null",breaks=50, xlab="ITT, percentage points")
abline(v=mean(0.01441675*100), col="red", lty=3)

sum((itt.store.sn)>(0.01441675))/sims


#FIG. AH4. RAND. DIST. OF ITT ON MEDIUM-HIGH-PROPENSITY VOTERS, SHARP NULL
#Estimate Model 7 without block fixed effects to generate point estimate for randomization inference
model7.noblock <- miceadds::lm.cluster( data = df.medhigh , formula = voucher ~ assignment, 
                                        cluster = "householdID" )
summary(model7.noblock)

set.seed(101)
sims <- 10000
itt.store.sn <- rep(NA,sims)

# Set stratification, balance variables, and creative fields
strat_vars <- c("votepropensity","multivoterhousehold","wacanmember", "momcat", "reachable_cellcat", "reachable_emailcat")
bal_vars <- c("agescale","sex","racename","incomescore","educationlevel")
#Subset to required fields
household.sub <- subset(df, select = c(strat_vars, bal_vars,"householdID","randomizationrecords","voucher","block"))
#Remove all but household representatives from the subset
household.sub <- household.sub[household.sub$randomizationrecords==1,]

# Compile randomization distribution of the ITT under the sharp null
for (i in 1:sims){
  household.sub$assign_treat<-NA
  household.sub$assign_treat<-block_ra(blocks=household.sub$block,prob=0.5)
  bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
  p<-(coef(summary(bal)))[,4]
  p<-as.numeric(p)
  while(p[1]<0.1|p[2]<0.1|p[3]<0.1|p[4]<0.1|p[5]<0.1|p[6]<0.1|p[7]<0.1|p[8]<0.1|p[9]<0.1|
        p[10]<0.1|p[11]<0.1|p[12]<0.1|p[13]<0.1|p[14]<0.1|p[15]<0.1|p[16]<0.1|p[17]<0.1|
        p[18]<0.1|p[19]<0.1|p[20]<0.1|p[21]<0.1|p[22]<0.1|p[23]<0.1){
    household.sub$assign_treat<-NA
    household.sub$assign_treat<-block_ra(blocks=household.sub$block)
    bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
    p<-(coef(summary(bal)))[,4]
    p<-as.numeric(p)
  }
  
  ##Link subset back to main DF
  #Create assignment variable, assign duplicates to the values of their designated housemate
  df$D <- NA
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==1], 1, df$D)
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==0], 0, df$D)
  #Drop households with 8 or more voters
  df<-df[!df$householdID==24341,]
  df<-df[!df$householdID==27505,]
  df<-df[!df$householdID==1341,]
  df<-df[!df$householdID==58,]
  df<-df[!df$householdID==10117,]
  df<-df[!df$householdID==15890,]
  df<-df[!df$householdID==17965,]
  df<-df[!df$householdID==13328,]
  df<-df[!df$householdID==3609,]
  df<-df[!df$householdID==28910,]
  df<-df[!df$householdID==371,]
  df<-df[!df$householdID==1327,]
  df<-df[!df$householdID==13751,]
  df<-df[!df$householdID==27364,]
  df<-df[!df$householdID==4832,]
  df<-df[!df$householdID==13641,]
  df<-df[!df$householdID==25362,]
  df<-df[!df$householdID==9485,]
  df<-df[!df$householdID==10331,]
  df<-df[!df$householdID==18711,]
  df<-df[!df$householdID==19278,]
  #Estimate ITT under sharp null iteratively, storing each estimate
  df.medhigh<-df[df$medhigh==1,]
  itt.store.sn[i] <- mean(df.medhigh$voucher[df.medhigh$D==1],na.rm=TRUE) - 
    mean(df.medhigh$voucher[df.medhigh$D==0],na.rm=TRUE)
}
hist(itt.store.sn*100, main="Fig. AH4. Rand. Dist. of ITT on Medium-High-Propensity Voters, Sharp Null",breaks=50, xlab="ITT, percentage points")
abline(v=mean(0.004176003*100), col="red", lty=3)

sum((itt.store.sn)>(0.004176003))/sims


#FIG. AH5. RAND. DIST. OF ITT ON MEDIUM-LOW-PROPENSITY VOTERS, SHARP NULL
#Estimate Model 6 without block fixed effects to generate point estimate for randomization inference
model6.noblock <- miceadds::lm.cluster( data = df.medlow , formula = voucher ~ assignment, 
                                        cluster = "householdID" )
summary(model6.noblock)

set.seed(101)
sims <- 10000
itt.store.sn <- rep(NA,sims)

# Set stratification, balance variables, and creative fields
strat_vars <- c("votepropensity","multivoterhousehold","wacanmember", "momcat", "reachable_cellcat", "reachable_emailcat")
bal_vars <- c("agescale","sex","racename","incomescore","educationlevel")
#Subset to required fields
household.sub <- subset(df, select = c(strat_vars, bal_vars,"householdID","randomizationrecords","voucher","block"))
#Remove all but household representatives from the subset
household.sub <- household.sub[household.sub$randomizationrecords==1,]

# Compile randomization distribution of the ITT under the sharp null
for (i in 1:sims){
  household.sub$assign_treat<-NA
  household.sub$assign_treat<-block_ra(blocks=household.sub$block,prob=0.5)
  bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
  p<-(coef(summary(bal)))[,4]
  p<-as.numeric(p)
  while(p[1]<0.1|p[2]<0.1|p[3]<0.1|p[4]<0.1|p[5]<0.1|p[6]<0.1|p[7]<0.1|p[8]<0.1|p[9]<0.1|
        p[10]<0.1|p[11]<0.1|p[12]<0.1|p[13]<0.1|p[14]<0.1|p[15]<0.1|p[16]<0.1|p[17]<0.1|
        p[18]<0.1|p[19]<0.1|p[20]<0.1|p[21]<0.1|p[22]<0.1|p[23]<0.1){
    household.sub$assign_treat<-NA
    household.sub$assign_treat<-block_ra(blocks=household.sub$block)
    bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
    p<-(coef(summary(bal)))[,4]
    p<-as.numeric(p)
  }
  
  ##Link subset back to main DF
  #Create assignment variable, assign duplicates to the values of their designated housemate
  df$D <- NA
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==1], 1, df$D)
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==0], 0, df$D)
  #Drop households with 8 or more voters
  df<-df[!df$householdID==24341,]
  df<-df[!df$householdID==27505,]
  df<-df[!df$householdID==1341,]
  df<-df[!df$householdID==58,]
  df<-df[!df$householdID==10117,]
  df<-df[!df$householdID==15890,]
  df<-df[!df$householdID==17965,]
  df<-df[!df$householdID==13328,]
  df<-df[!df$householdID==3609,]
  df<-df[!df$householdID==28910,]
  df<-df[!df$householdID==371,]
  df<-df[!df$householdID==1327,]
  df<-df[!df$householdID==13751,]
  df<-df[!df$householdID==27364,]
  df<-df[!df$householdID==4832,]
  df<-df[!df$householdID==13641,]
  df<-df[!df$householdID==25362,]
  df<-df[!df$householdID==9485,]
  df<-df[!df$householdID==10331,]
  df<-df[!df$householdID==18711,]
  df<-df[!df$householdID==19278,]
  #Estimate ITT under sharp null iteratively, storing each estimate
  df.medlow<-df[df$medlow==1,]
  itt.store.sn[i] <- mean(df.medlow$voucher[df.medlow$D==1],na.rm=TRUE) - 
    mean(df.medlow$voucher[df.medlow$D==0],na.rm=TRUE)
}
hist(itt.store.sn*100, main="Fig. AH5. Rand. Dist. of ITT on Medium-Low-Propensity Voters, Sharp Null",breaks=50, xlab="ITT, percentage points")
abline(v=mean(0.003376291*100), col="red", lty=3)

sum((itt.store.sn)>(0.003376291))/sims


#FIG. AH6. RAND. DIST. OF ITT ON LOW-PROPENSITY VOTERS, SHARP NULL
#Estimate Model 5 without block fixed effects to generate point estimate for randomization inference
model5.noblock <- miceadds::lm.cluster( data = df.lowprop , formula = voucher ~ assignment, 
                                        cluster = "householdID" )
summary(model5.noblock)

set.seed(101)
sims <- 10000
itt.store.sn <- rep(NA,sims)

# Set stratification, balance variables, and creative fields
strat_vars <- c("votepropensity","multivoterhousehold","wacanmember", "momcat", "reachable_cellcat", "reachable_emailcat")
bal_vars <- c("agescale","sex","racename","incomescore","educationlevel")
#Subset to required fields
household.sub <- subset(df, select = c(strat_vars, bal_vars,"householdID","randomizationrecords","voucher","block"))
#Remove all but household representatives from the subset
household.sub <- household.sub[household.sub$randomizationrecords==1,]

# Compile randomization distribution of the ITT under the sharp null
for (i in 1:sims){
  household.sub$assign_treat<-NA
  household.sub$assign_treat<-block_ra(blocks=household.sub$block,prob=0.5)
  bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
  p<-(coef(summary(bal)))[,4]
  p<-as.numeric(p)
  while(p[1]<0.1|p[2]<0.1|p[3]<0.1|p[4]<0.1|p[5]<0.1|p[6]<0.1|p[7]<0.1|p[8]<0.1|p[9]<0.1|
        p[10]<0.1|p[11]<0.1|p[12]<0.1|p[13]<0.1|p[14]<0.1|p[15]<0.1|p[16]<0.1|p[17]<0.1|
        p[18]<0.1|p[19]<0.1|p[20]<0.1|p[21]<0.1|p[22]<0.1|p[23]<0.1){
    household.sub$assign_treat<-NA
    household.sub$assign_treat<-block_ra(blocks=household.sub$block)
    bal<-glm(assign_treat~agescale+sex+racename+incomescore+educationlevel, data=household.sub, family=binomial(link=logit))
    p<-(coef(summary(bal)))[,4]
    p<-as.numeric(p)
  }
  
  ##Link subset back to main DF
  #Create assignment variable, assign duplicates to the values of their designated housemate
  df$D <- NA
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==1], 1, df$D)
  df$D <- ifelse(df$householdID %in% household.sub$householdID[household.sub$assign_treat==0], 0, df$D)
  #Drop households with 8 or more voters
  df<-df[!df$householdID==24341,]
  df<-df[!df$householdID==27505,]
  df<-df[!df$householdID==1341,]
  df<-df[!df$householdID==58,]
  df<-df[!df$householdID==10117,]
  df<-df[!df$householdID==15890,]
  df<-df[!df$householdID==17965,]
  df<-df[!df$householdID==13328,]
  df<-df[!df$householdID==3609,]
  df<-df[!df$householdID==28910,]
  df<-df[!df$householdID==371,]
  df<-df[!df$householdID==1327,]
  df<-df[!df$householdID==13751,]
  df<-df[!df$householdID==27364,]
  df<-df[!df$householdID==4832,]
  df<-df[!df$householdID==13641,]
  df<-df[!df$householdID==25362,]
  df<-df[!df$householdID==9485,]
  df<-df[!df$householdID==10331,]
  df<-df[!df$householdID==18711,]
  df<-df[!df$householdID==19278,]
  #Estimate ITT under sharp null iteratively, storing each estimate
  df.lowprop<-df[df$lowprop==1,]
  itt.store.sn[i] <- mean(df.lowprop$voucher[df.lowprop$D==1],na.rm=TRUE) - 
    mean(df.lowprop$voucher[df.lowprop$D==0],na.rm=TRUE)
}
hist(itt.store.sn*100, main="Fig. AH6. Rand. Dist. of ITT on Low-Propensity Voters, Sharp Null",breaks=50, xlab="ITT, percentage points")
abline(v=mean(0.001265398*100), col="red", lty=3)

sum((itt.store.sn)>(0.001265398))/sims

#TABLE AI1. HETEROGENEOUS TREATMENT EFFECTS ON VOUCHER USE BY VOTING HISTORY
#Treatment effect on low-propensity voters
#Count observations for Model 5
model.lpm.lowprop<-lm(voucher~assignment+block,data=df.lowprop)
nobs(model.lpm.lowprop)
#Model 5
model5 <- miceadds::lm.cluster( data = df.lowprop , formula = voucher ~ assignment+block, 
                                cluster = "householdID" )
summary(model5)

#Treatment effect on medium-low-propensity voters
#Count observations for Model 6
model.lpm.medlow<-lm(voucher ~ assignment+block,data=df.medlow)
nobs(model.lpm.medlow)
#Model 6
model6 <- miceadds::lm.cluster( data = df.medlow , formula = voucher ~ assignment+block, 
                                cluster = "householdID" )
summary(model6)

#Treatment effect on medium-high-propensity voters
#Count observations for Model 7
model.lpm.medhigh<-lm(voucher~assignment+block,data=df.medhigh)
nobs(model.lpm.medhigh)
#Model 7
model7 <- miceadds::lm.cluster( data = df.medhigh , formula = voucher ~ assignment+block, 
                                cluster = "householdID" )
summary(model7)

#Treatment effect on high-propensity voters
#Count observations for Model 8
model.lpm.highprop<-lm(voucher~assignment+block,data=df.highprop)
nobs(model.lpm.highprop)
#Model 8
model8 <- miceadds::lm.cluster( data = df.highprop , formula = voucher ~ assignment+block, 
                                cluster = "householdID" )
summary(model8)

#TABLE AJ1. ITT EFFECTS OF VOTER MOBILIZATION ON VOTER TURNOUT
#Count observations for Model 9
model.lpm.voted<-lm(voted~assignment+block, data=df)
nobs(model.lpm.voted)
#Model 9
model9 <- miceadds::lm.cluster( data = df , formula = voted ~ assignment+block, 
                                cluster = "householdID" )
summary(model9)
#Count observations for Model 10
model.lpm.voted.ctrl<-lm(voted ~ assignment+educationscale+incomescale+age+female+black+asian+latinx+block,data=df)
nobs(model.lpm.voted.ctrl)
#Model 10
model10 <- miceadds::lm.cluster( data = df , formula = voted ~ assignment+educationscale+incomescale+age+female+black+asian+latinx+block, 
                                 cluster = "householdID" )
summary(model10)
#Count observations for Model 11
model.lpm.voted.int<-lm(voted ~ assignment+member+highprop+assignment*member+assignment*highprop+block,data=df)
nobs(model.lpm.voted.int)
#Model 11
model11 <- miceadds::lm.cluster( data = df , formula = voted ~ assignment+member+highprop+assignment*member+assignment*highprop+block, 
                                 cluster = "householdID" )
summary(model11)
#Count observations for Model 12
model.lpm.voted.int.ctrl<-lm(voted ~ assignment+member+highprop+assignment*member+assignment*highprop+educationscale+incomescale+age+female+black+asian+latinx+block,data=df)
nobs(model.lpm.voted.int.ctrl)
#Model 12
model12 <- miceadds::lm.cluster( data = df , formula = voted ~ assignment+member+highprop+assignment*member+assignment*highprop+educationscale+incomescale+age+female+black+asian+latinx+block, 
                                 cluster = "householdID" )
summary(model12)


#TABLE AK1. HETEROGENEOUS TREATMENT EFFECTS ON VOTER TURNOUT BY VOTING HISTORY
#Count observations for Model 13
model.lpm.vote.lowprop<-lm(voted~assignment+block,data=df.lowprop)
nobs(model.lpm.vote.lowprop)
#Model 13
model13 <- miceadds::lm.cluster( data = df.lowprop , formula = voted ~ assignment+block, 
                                 cluster = "householdID" )
summary(model13)
#Count observations for Model 14
model.lpm.vote.medlow<-lm(voted~assignment+block,data=df.medlow)
nobs(model.lpm.vote.medlow)
#Model 14
model14 <- miceadds::lm.cluster( data = df.medlow , formula = voted ~ assignment+block, 
                                 cluster = "householdID" )
summary(model14)
#Count observations for Model 15
model.lpm.vote.medhigh<-lm(voted~assignment+block,data=df.medhigh)
nobs(model.lpm.vote.medhigh)
#Model 15
model15 <- miceadds::lm.cluster( data = df.medhigh , formula = voted ~ assignment+block, 
                                 cluster = "householdID" )
summary(model15)
#Count observations for Model 16
model.lpm.vote.highprop<-lm(voted~assignment+block,data=df.highprop)
nobs(model.lpm.vote.highprop)
#Model 16
model16 <- miceadds::lm.cluster( data = df.highprop , formula = voted ~ assignment+block, 
                                 cluster = "householdID" )
summary(model16)
